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The recently reported ultracold trapped Fermi gases with tunable atomic scattering 
length [1-11] in the vicinity of a Feshbach resonance stimulated a large number of theo- 
retical investigations. Some of these works are based on the assumption that the properties 
of strongly interacting dilute Fermi gas at zero temperature are well described by the hydro- 
dynamic approximation (HA) [12-15] 

^ + V(nv) = 0, (1) 

| + I V (K„ + ^)) + i^)=0, (2) 
at m on 2 

where n is the density, e(n) is the ground-state energy per particle of the homogeneous 
system and v is the velocity field. 

In this paper the dynamics of strongly interacting trapped dilute Fermi gases (dilute in 
the sense that the range of interatomic potential is small compared with inter-particle spacing 
) is investigated in the single equation approach to the time-dependent density-functional 
theory covering the whole crossover region at zero temperature. It is shown for the case of 
elongated cigar-shaped harmonic traps that the calculated corrections to the HA may be 
important even for systems with a rather large number of atoms. 

We mention here Refs.[16] where an extension of the density-functional theory (DFT) to 
superconducting systems [17] was generalized to a number of nuclear and atomic systems. 

Let us consider a Fermi gas consisting of a 50-50 mixture of two different states con- 
fined in a harmonic trap V ext (r) = (m/2)(uj_(x 2 + y 2 ) + uj 2 z 2 ). In Eq.(2), the kinetic- 
energy density t(n) is approximated by the Thomas-Fermi (TF) kinetic-energy density 
^tf(^) = (3/l0)nh 2 kp/m, where hp = (S^n) 1 ^ 3 . For slowly varying densities character- 
ized by the condition | Vn | /n 4 / 3 <C 1, the kinetic energy density is well represented 
by the Kirzhnitz gradient expansion (KGE) [18] t(n) = txF{n) + tw{n)/9 + where 
tw( n ) — (^ 2 /(8m))(Vn) 2 /n is the original von Weizsacker density (OWD)[19], which gives 
the entire kinetic energy density of noninteracting bosons. 

In the case of large but finite number of atoms N, the density n is not constant. At small 
distances the ratio | Vn | /n 4 / 3 is small and both the Kirzhnitz correction and the OWD 
are negligible. On the contrary, near the surface the Hartree-Fock (HF) type densities are 
proportional to the square of the last occupied state. Therefore, the OWD is important in 
this case and it is expected to determine the asymptotic behavior of the density at large 



2 



distances. It is also expected that the OWD is important in the case of the tight radial 
trapping, A< 1. In Refs.[20], the OWD was considered as a correction to the TF kinetic- 
energy density. 

Adding the OWD to £tf(^) we have 

dv l„, Tr d(ne(n)) 1 9 h 2 1 ~ 

m + ^ (v '« + + 2 mv ~ 2^ vV5) = a (3) 

We define the density of the system as n(f,t) =\ *&(r,t) | 2 , and the velocity field v as 
v(f,t) = 7i(#*Vtf - tfVtf*)/(2imn(f,f)). From Eqs.(l) and (3), we obtain the following 
nonlinear Schrodinger equation 

._dV h 2 2T d(ne(n)) T 

which is equivalent, to a certain extent, to the single equation approach of Deb et al. [21] to 
the time-dependent density- functional theory (TDDFT). 

If the trap potential, V ex t, is independent of time, one can write \l/(r, t) = $(f) exp(—if/t/h), 
where \x is the chemical potential, and $ is normalized to the total number of particles, 
J df\ $ | 2 = N. Then Eq.(4) becomes 

/ r-,i t t 9(ne(n)) .^ _ . . 

where the solution of the equation (5) minimizes the energy functional E = N < $ | 
— |^V 2 + V ext + e(n) | $ >, and the chemical potential \x is given by /i — dE/dN. 

In order to take into account atoms lost by inelastic collisions, we model the loss by the 
rate equation 

where xi^jt) = J2i=i h^Qiij 1 ), n l gi is the local /-particle correlation function and k\ is the 
rate constant for the /-body atoms loss. The generalization of Eq.(4) for the case of inelastic 
collisions reads [22] 

..<9\l/ ^ 2 „9 T Tr T d(ne(n)) T .^v^, / i , \ T 

For the negative S-wave scattering length between the two fermionic species, a < 0, in 
the low-density regime, fcf | o |< 1, the ground state energy per particle , e(n), is well 
represented by an expansion in power of kp \ a \ [26] 

3 1 

e(n) = 2E F \— - —k F I a I +0.055661(A; F I a I) 2 - 0.00914(A; F I a I) 3 + ...], (7) 
10 Jff 



where E F = h 2 k F /(2m). In the opposite regime, a — > — oo (the Bertsch many-body problem, 
quoted in Refs.[27]), e(n) is proportional to that of the non-interacting Fermi gas 

where a universal parameter f3 is estimated to be (3 = —0.56 [28]. 

In the a — > +0 limit the system reduces to the dilute Bose gas of dimers [29] 

e(n) = E F {-l/{k F a) 2 + a m k F /(QTx) + ) (9) 

where a m is the boson-boson scattering length. While the BCS mean-field theories [30] 
predict a m = 2a [31], a solution of 4-fermion problem for contact scattering provided the 
value a m ~ 0.6a [32]. 

Very little is known about the correct form of e(n) in the intermediate range. Therefore, 
a simple interpolation of the form e(n) m E F P{k F a) with a smooth function P(x) mediating 
between the known limits suggests itself as a pragmatic alternative. 

In Ref.[33] it has been proposed a [2/2] Pade approximant for the function P(x) for the 
negative a 

p(,) = ^- 2 / 1 ' 3 ; l ^f 2 , (io) 
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where Si = 0.106 1 03 , 5 2 = 0.1875 1 5 , 5 3 = 2.29 188 , 5 4 = 1.11616. Eq.(10) is constructed 
to reproduce the first four terms of the expansion (6) in the low-density regime and also 
to reproduce exactly results of the recent Monte Carlo calculations [28], (3 = —0.56, in the 
unitary limit, k F a — > — oo. 

For the positive a case ( the interaction is strong enough to form bound molecules with 
energy E mo i) we consider a [2/2] Pade approximant 

p( r \ = Emo1 i aiX + 012x2 h-n 

{) 2E F + l + a 3 x + a 4 x 2 ' 1 J 

where parameters a are fixed by two continuity conditions at large x, 1/x — > 0, and by 
two continuity conditions at small x. For example, a.\ = 0.0316621, a 2 = 0.0111816, 0:3 = 
0.200149, and a A = 0.0423545 for a m = 0.6a. 

Fig. 1 and Fig. 2 show the comparison between [2/2] Pade approximations, Eqs.(10,ll), 
and the lowest order constrained variational (LOCV) approximation [34] and the BCS mean- 
field theory for e(n). The LOCV calculations agree very well with the [2/2] Pade approxi- 
mation results on the BCS side (a < 0). It is evident the difference between our results and 



the BCS mean-field theory calculations. For example, the BCS mean-field gives (3 = —0.41. 
We mention here that e(n)/E F on the BCS side (a < 0) and (e(n) + \E mol \/2) / E F on the 
BEC side (a > 0) show a smooth monotonic behavior as a function of k F a. 

The predictions of Eq.(5) with e(n) from Eq.(10) for the axial cloud size of strongly 
interacting 6 Li atoms are shown in Fig 3 [35]. It indicates that the TF approximation of 
the kinetic energy density is a very good approximation for the experimental conditions of 
Ref.[ll], NX « 10 4 (inclusion of the OWD gives a negligible effect, < 0.5%) [37]. 

It can be proved [24] that every solution of equation (4) is a stationary point of an action 
corresponding to the Lagrangian density 

C = - (*— - * *_) + _ I W | 2 + e(n)n + V ext n, 

which for \1/ = e l ^ r ' t ^n 1 ^ 2 (f,t) can be rewritten as 

h 2 h 2 
C = Hn + ^(Vv 7 ^) 2 + ^ n ( V< ^) 2 + e ( n ) n + V ^ n - 

For a time-dependent harmonic trap, V ext (f, t) = (m/2) X)f=i ^fi^x 2 , a suitable trial function 
can be taken as (f>(r,t) = x(t) + {m/{2fi) Yh=i Vi{t) x l > n(f,i) = rio(xi/bi(t))/((t), where 
C(0 = Ilj With this ansatz, the Hamilton principle, 5 J dt J C d 3 r = 0, gives the following 
equations for the scaling parameters bi 

2 < T > 1 r 

h + u 2 ^ - - - B - m<x 2 >bi J [n 2 de(n)/dn] n=noim{t) d 3 r((t) = 0, (12) 

where 6i(0) = 1, 6«(0) = and a>j = u>i(0) fix the initial configuration of the system, cor- 
responding to the density no(f) and < Tj >= —% 2 /{2mN) J n}l 2 {d 2 /dx^r^^dPr •, < x 2 >= 
(1/iV)/ nx 2 ci 3 r. 

Expanding Eqs.(12) around equilibrium (6j = 1) we get the following equations for the 
collective frequencies, u> 

uo 2 1 

(2 + m- —)vi + (i + -Ki + Xi)(yi + V2 + ys) = o, (13) 

where Ki = 4 < > /(muo 2 < x 2 >) and Xi — S ^o<9 2 e/ ( y dn Q l )d 3 r / (muo 2 < x 2 >). 

In Table I we give the calculated values of the radial breathing mode frequency, v = 
^rad/(27r), of highly degenerate gas of 6 Li atoms near a Feshbach resonance at 822 G [41]. 
It can be seen from Table I, that the difference between two approximations, Eqs.(l,2) 
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and Eq.(4), is less than 0.75%, and both approximations give a very good agreement with 
experimental data of Ref.fl]. The parameter AiV for this case is very large, AiV > 10 4 . 

In Fig. 4, we present the calculations for the frequency of the radial compression mode 
oj r ad as a function of the dimensional parameter (A^ 1//6 a/a/ lo ) _1 in the case of an anisotropic 
trap (u x = uj y = w±, uj z /uj± = A). One can easily see that the corrections to the hydrody- 
namic approximation (HA), Eqs.(l) and (2), are important even for relatively large N and 
XN. For example, the correction to u ra d in unitary limit is larger than 11% and 25% for 
A = 1(T 2 , N = 10 4 and A = 1(T 2 , N = 10 3 , respectively. 

In the HA, oo rad is independent of N for a fixed (N l l 6 a/ a ho )~ l . The deviation from this 
behavior does not demonstrate the cross-over to the ID behavior, since AiV > 1 [42]. It 
demonstrates that the validity of the HA depends on the properties of the trap. In Ref.[43] 
it was shown that, for the case of isotropic trap, A = 1, with N = 20 and N = 240, the TF 
approximation reproduces the energy within accuracies of 2% and 1%, respectively. 

In Fig. 5, the calculated radial compressional frequency is compared with experimental 
data [1] in the BCS-BEC crossover region. There is a very good agreement between calcula- 
tions and experimental data [1]. However our calculations for uw disagree with experimental 
data of Ref.[44]. 

In the present paper, we have used Eq.(4). The next step is to develop the Kohn-Sham 
time-dependent DFT [45] for two-component Fermi gases in elongated traps (A < 1), which 
we will consider in our future work. 

In conclusion, the dynamics of strongly interacting trapped dilute Fermi gases is inves- 
tigated in the single equation approach to the time-dependent density-functional theory. 
Our results are in good agreement with recent experimental data in the BCS-BEC crossover 
regime. 

It is also shown that the calculated corrections to the hydrodynamic approximation may 
be important even for systems with a rather large number of atoms. 
ALZ thanks V. Dunjko for useful discussion. 

Note added. - While this work was being prepared for publication, preprints [46] appeared 
in which the authors calculate the equation of state, e(n), using the quantum Monte Carlo 
method. Their results are in a good agreement with ours Pade [2/2] approximation for both 
negative and positive scattering length. 
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Table I. Radial breathing mode frequencies v = uj rad /(2-K) of highly degenerate gas of 
6 Li atoms near a 822(G) Feshbach resonance [36]. B is applied magnetic field, u exp indicate 
experimental data from the Duke University group [1], v and vtf represent theoretical 
calculations that use equation (4) and the hydrodynamic approximation (equations (1) and 
(2)), respectively. The trap parameters are u± = 2n x 1549, uj z = 2n x 70. 
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Fig.l. The ground state energy per particle, e(n), in units of 3h 2 k 2 F /(10m) function 
of the gas parameter (fcpa) -1 . The solid line, the long dashed line and the short dashed 
line represent the results calculated using the [2/2] Pade approximation, Eq.(10), the LOCV 
approximation, and the BCS mean-field theory, respectively. 
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Fig. 2. The ground state energy per particle, e(n) + \E mol \/2, in units of Ti 2 k 2 F j '(2m) as a 
function of the gas parameter (kpa)^ 1 . The dashed line, the dotted-dashed line and the 
solid line represent the results calculated using the BCS mean-field theory, the [2/2] Pade 
approximation, Eq.(ll), with a m = 2a, and a m = 0.6a, respectively. 
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Fig. 3. Axial cloud size of strongly interacting 6 Li atoms after normalization to a non- 
interacting Fermi gas with N = 4 x 10 5 atoms as a function of the magnetic field B [32]. The 
trap parameters are cu± = 2tc x 640Hz, oj z = 2n(600B /kG + 32) 1 / 2 Hz. The solid line and 
dashed line represent the results of theoretical calculation that includes the OWD or uses 
the TF approximation for the kinetic energy density, respectively. The circular dots indicate 
experimental data from the Innsbruck group [11]. 
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Fig. 4. Radial compressional frequency, uw, of the cloud of the iV = 10 4 fermions (solid 
line) and iV = 10 3 fermions (dashed line) in unit of uj± as a function of the dimensional 
parameter (iV 1 / 6 a/a/ lo ) _1 . The trap parameter A is assumed to be equal to 10~ 2 . The lower 
line (dashed-dotted line) represents the results in the hydrodynamic approximation, Eqs. (1) 
and (2), in which cu rac i is independent of N for a fixed (iV^a/a^) -1 . 
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Fig. 5. The radial compressional frequency as a function of (A^ 1 / 6 a/a/ lo ) _1 . The solid line 
and the dashed line represent the results calculated using the [2/2] Pade approximation with 
a m = 0.6a and a m = 2a, respectively. The solid circles with error bars are the experimental 
results given by the Duke University group [1]. 
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